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I . INTRODUCTION 



The main concern when approximating either an ordinary or a partial 
differential equation by a discrete equation, either finite difference or 
finite element, is the degree to which the discrete equation solution 
approximates the solution to the differential equation. The closeness of 
this approximation has traditionally been considered from both its 
quantitative aspects, such as relative error (Henrici (1963)), and its 
qualitative aspects, such as behavior of transients, propagation of fronts, 
etc. (Haltiner (1971)). Analyses based primarily on consideration of the 
basically qualitative property of numerical dispersion, where the 
discretization process causes distortions in the propagation velocities 
for different frequencies, are especially prevalent (Arakawa and Lamb 
(1977)). 

In this paper, we apply certain ideas derived from the concepts of 
transfer functions and digital filters in electrical engineering to the 
study of the qualitative behavior of discretization schemes. The 
particular model we shall use is based on the linearized, one dimensional 
shallow water equations, without mean flow, a model which has important 
application to the study of the process of geostrophic adjustment. (It 
will be evident that the basic approach, however, is not limited to this 
model.) The qualitative effects of discretization schemes on this 
equation have been studied by Winninghoff (1968), Arakawa and Lamb (1977), 
Schoenstadt (1977), and others. We shall show that use of the transfer 
function concept leads to important insights into the differences caused 
by different discretization schemes - differences that are not fully 
evident from phase propagation considerations only. 
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II. TRANSFER FUNCTIONS 

It is well known (Stremler (1977)) that linear, space (or time) 
invariant systems can be fully described in terms of their so-called 
transfer function. (Figure 1) That is, if we denote the input to the 
system as y^(x), and the output as y Q (x), then the Fourier transforms of 
the input and output, denoted respectively as 



a, 

y^k) 



J r y^(x) e '*'^ CX dx 







(k) 



J y Q (x) e dx 



are related by 



( 1 ) 



y o (k) = ?(k)y.(k) . (2) 

Equation (2) is commonly said to represent the system in the transform 
domain. The equivalent representation in the physical (time or space, as 
appropriate) domain is 

f 00 

y Q (x) = / <K X - s) y ± (s) ds . (3) 

In equation (3), (f>(x) is commonly referred to as the impulse response of 
the system, and is interpreted as the response of the system to an input 
delta function at x = 0, i.e. 6(x). 

In general, 4> (k) is a complex valued function. The physical interpre- 

% 

tation of both the magnitude and phase of cf> (k) is easily arrived at by 
considering the response of the system to a single sinusoidal input of 
arbitrary frequency, i.e. 

y.(x) = e lkiX . (4) 
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It is very straightforward to show that 

ik .x 

y Q (x) = ^(k ± ) e 1 

Thus, if we represent 4>(k) in polar form as 



?(k.) 




equation (5) becomes: 



> 



i(k.x + k.) 

y Q (x) = |^(k ± ) | e 



ik (x + (k ,/k ) ) 
- I?(k 1 )| e 1 6 1 



(5) 



Clearly, from (5), we can interpret the magnitude of the transfer function 
as the factor by which the amplitude of a sinusoid of a given frequency is 
either amplified or attenuated, and the argument of the transfer function 
as determining the amount by which the phase of the output sinusoid is 
shifted relative to the phase of the input sinusoid. 

In general, neither |(f>(k) | nor (k^/k) is constant, and hence, the 
different frequencies in a given input are amplified/attenuated and shifted 
in phase by different amounts. Thus the output signal generally has its 
shape (graph) altered from the input signal. (Figure 2) This alteration 
of shape is commonly called distortion, and, based on our discussion, is 
composed of the effects of two actions - commonly called amplitude 
distortion, which is due to the deviation of |<{>(k) | from a constant, and 
phase (or delay) distortion, due to the deviation of k^ from a linear 
function of k. Clearly, to completely understand the effect of a linear 
system on an input signal, one must know both of these effects. 
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In electrical engineering, it is common to call any linear, time 
(or space) invariant device, which can be described by a transfer function 
as discussed above, a filter. We shall use the terminology hereafter in 
this paper. 

Lastly we note that the energy in an "output" is proportional to 



and thus one can also use the transfer function to describe the relation 
between the energies contained in the "input" and "output" to a linear. 





( 6 ) 





dx 



• * 2 

invariant system. (Note the quantity |y(k)| is commonly referred to as 



the spectral density of the signal y(x).) 
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III. DISCRETE FILTERS AND TRANSFER FUNCTIONS 



A great deal of research is presently directed toward the study of 
so-called discrete, or digital filters, which arise when the continuous 
processes, such as differentiation and integration, in an invariant linear 
filter are replaced by analogous approximations based on sampled, or 
discrete data. Thus, for example, 

y(x + Ax ) = y (x) + f (x) (Ax) , (7) 

is a digital filter approximating, in some sense (specifically by using 
an Euler forward predictor), the continuous integrator: 

^ = f(x) • (8) 
The transfer functions of such filters can be easily determined using 
Fourier transforms, and the observation that 

y(x 4- Ax ) e dx = e y(k) . (9) 




Thus, the transfer function of the process given by (7) is (where f(x) is 
considered the f, input n and y(x) the "output") , after some algebra, 



$(k) = 



(ax)e ik(ax)/2 

2i sin(k(Ax)/2) 



(Note the transfer function of the continuous process is <f> (k) = — .) 

There is, however, one significant difference between the transfer 
functions of a continuous filter and its digital (discrete) analog. This 
difference is commonly referred to with the terms aliasing or folding, and 
the basic work on it is commonly attributed to Nyquist (Clark (1977)). In 
its simplest interpretation, Nyquist T s work implies that a sampling device, 
sampling at intervals of ( Ax ) , is incapable of resolving waves of 
frequency greater than (tt/(Ax)), and, if any energy is actually present 
in a sampled continuous signal at these higher frequencies, it will be 
erroneously resolved (aliased) into a frequency lower than (tt/(Ax ))• 
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Thus, while a continuous filter can, in principle, react to any input 
frequency, and hence has a transfer function defined for all k, that of 
a digital filter is considered only for (complex) frequencies between 
+(tt/(Ax)), the so-called Nyquist cut-off or Nyquist limit. 

Lastly, we note that if a digital filter has a significant amplitude 
response (i.e. if |$(k) | is not close to zero) near the Nyquist limit, 
then the output of that system may have a significant portion of its 
energy in these high frequency oscillations. These oscillations often 
appear to the observer as noise, or u noise-like !! , and hence, in discussing 
digital filter design, many books and articles recommend controlling the 
transfer function response near this limit. (This is sometimes referred 
to as "windowing 11 . ) 
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IV. THE SHALLOW WATER EQUATIONS AS A FILTER 



The basic equations we propose to study, using the ideas discussed 
above, are the linearized, one-dimensional, shallow water equations in 
an infinite region, with no mean flow: 



3u 3h 

' fv + 8 " 0 



3v 

at 

ah 



+ fu = 0 
3u 



( 10 ) 



3t + 1 *S-° ’ 



where u denotes the perturbation velocity in the x direction, v the 
perturbation velocity normal to the x direction, H and h the mean and 
perturbed heights of the free surface, respectively, and g > 0 and f > 0 
the gravitational and Coriolis parameters, respectively. This model is 
especially important in the study of the meteorological process called geo- 
strophic adjustment, which has been studied in some detail from several 
approaches by Rossby (1937), Cahn (1945), Winninghoff (1968), Blumen (1972), 
Arakawa and Lamb (1977) and Schoenstadt (1977) . This process of geostrophic 
adjustment is important because it is the primary mechanism by which an 
atmosphere, modeled by the so-called meteorological primitive equations, 
reacts to errors (either numerical or due to errors in observational data) 
in the initial data. The dispersive wave nature of (10) is the primary 
mechanism by which these errors (commonly referred to as "imbalances") event- 
ually spread out over the domain of interest, until the system reaches a so- 
called "balanced" or "geostrophically adjusted" condition. Forecasts at times 
before this balance is reached are generally inaccurate. Thus a primary con- 
sideration in numerically modeling (10) is to insure that the numerical 
balance mechanism both fairly accurately models the physical balancing pro- 
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cedure, and that it proceeds relatively quickly. 

As discussed in Schoenstadt (1977) , equations (10) can be easily 
solved by using spatial Fourier transform. If we continue to denote 
Fourier transforms by an overhead tilde, e.g.: 



u(k,t) = / u(x,t)e dx , 



(ID 



etc. Then (10) reduces directly to: 



du 



— = fv - ikg h , 



dv 



dt 

dh 



dt 



= -f u , 

= -ik H u , 



( 12 ) 



together with initial conditions: 



u^ = u(k,0) = / u(x,0)e dx , 



(13) 



etc. Then since (12) represents simply a set of coupled, constant coefficient 
(in the sense of t) ordinary differential equations, their solution can be 
easily shown, by the usual eigenvalue/eigenvector approach, to be: 



f v . ikg h . 

u(k,t) = u. cos vt H sin vt sin vt , 

1 v v 



v(k,t) u. sin vt + 

vi 



lk 2 gH , f 2 ~ ikgf 

, 2 * ~2 COS Vt ( V i + 2 

V V J V 



jl - cos vtj In , 



(14) 



\ ikH ~ . _ ikHf 

n(k,t) u. sin vt - — r— 

vi 2 



/f 2 


2 ) 


< 1 - COS Vt > V. + < — 


+ ^cos vt 


l f i ) 2 


2 l 


v 1 ( v 





where : 



v 2 = f 2 + k 2 gH = f 2 (l + A 2 k 2 ), \ = /iH/f 



(15) 



We now consider (14) in view of our discussion of filters and transfer func- 
tions. Observe that if we denote 



y i = 



h. 

1 



and 



y o = 



(16) 
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respectively, as the "input" and "output" to the "system" described by 
(10), then we can express (14) as 



y Q = <j)(k > t)y i 



(17) 



where, 



<f>(k,t) = 



cos vt 



- — sin vt 

v 



v 



sin vt 



k 2 gH f 2 

— 2 ~ + ~2 cos vt 
v v 






sin vt 



^ 2 ^ ^ ” cos 

v 



ikH 

sm vt 



ikHf 



{1 - cos vt} < —r- 



f 2 k 2 gH I 

— 2 + — 2 cos vt 

v v ! 



Immediate comparison of (17)-(18) with (2) is not possible since (k , t ) 

depends not only on the transform variable ("k"), but also on one of the 

physical variables ("t"). However, since sin vt and cos vt can be 

+ iv t 

expressed in terms of e , and since k is real, it is easily seen 

that the sinusoidal terms in (18) express only the phase relationships 
(which clearly change with time), while the other, time independent, 
coefficients reflect the amplitude effects of the "filter". In fact, it is 
fairly easy to show (Schoenstadt (1977)) that the sinusoidal terms in (18) 
produce dispersive, transient waves, propagating throughout the region. 

(These waves are clearly dispersive since (v(k)/k) is not constant.) The 
time independent coefficients in (18) still function as amplitude distortion 
terms in the sense that they redistribute the energy that is in the input 
(initial condition) waves, in a manner that varies with wave number, into the 
output (solution) waves. 

A close examination of (18) shows that, in fact, the "amplitude distor- 
tion" in this "system" is really governed in each term by precisely one of 



U8) 



> 9 - 



the three factors. 



1 , 



v 



k 

v 



or 



k_ 

2 



9 



V 

or the square of one of these. The amplitude response curves for these three 
expressions are shown in Figure 3. It is clear that the terms with the co- 
efficient (1/v) have their low frequencies least affected, the terms with 

coefficient (k/v) their high frequencies least affected, and those with 
2 

coefficient (k/v ) their frequencies in some middle range least effected. 
Such responses are often referred, respectively, to as low pass, high pass, 
and band pass filters. 

Since the "system" described by (10) involves both spatial and temporal 
variation, it is well known that not only the phase velocity (v(k)/k) , but 
also the group velocity (t|~) is important in describing the behavior, since 
it is the latter which governs the rate at which energy actually propagates. 
Thus, in Figure 4, we display both of these velocities for the system (10). 

We note that additional insight into some of the qualitative behavior 
of the solution to (10) can be obtained by decomposing (18) into its transient 
and steady state components, i.e., writing 



/v 




*3 y± + <J> T y± 






/v 



(19) 



where <j> and <}> can be obtained by inspection, since only time dependent 
^ 1 

terms contribute to the transients, however, we shall not pursue this analysis 
here. 



As a final comment, note that the total energy of the system represented 
by (10) is given by: 



2 

u 



+ V 




Thus the spectral density of the steady state field can be shown to 





v 




be 



( 20 ) 
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Since the term (k/v) acts as a high-pass filter, and the term (f/v) as a 
low pass-filter, we see that the energy in the steady state field is deter- 
mined by a combination depending primarily on the high frequency energy in 
the initial v field, and the low frequency energy in the initial h field. 

We shall not proceed any further with this analysis, since our main 
interest is to analyze how discrete schemes, and the corresponding "filters" 
they produce, approximate the continuous (differential) model. However, 
there are certainly other insights into the general behavior of the physical 
model still possible from this filter point of view. 
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V. TRANSFER FUNCTIONS FOR QUASI- DISCRETE FINITE DIFFERENCE AND FINITE 
ELEMENT METHODS 

In this section, we extend the analysis of our continuous model to 
the quasi-discrete case, where continuous spatial operations are replaced 
by operations only at sample points, while the continuous time operations 
are retained. This approach follows that in Winninghoff (1968) and Arakawa 
and Lamb (1977). However, we extend our analysis to include a more general 
case of quasi-discrete schemes than they considered. Specifically, we 
define an invariant linear discretization operator about the point x as 
any operator having the form 

n 2 

L [ f (x) ] = l c f(x + n(Ax) ) , (21) 

n 

n=-n^ 

where the c denote real constants, 
n 

The discretization schemes we shall investigate are specifically those 
satisfying the following conditions: 

(a) All discretization operators are linear and time and space invariant. 

(b) All discretization operators operate only on symmetric sets of 
sample points, i.e. in (21), n^ = n^ • 

(c) The discretization operators are balanced, i.e. in (21) | c^ | = l c _ n l • 

(d) The operators are applied to (8)-(10) consistently, i.e. if one of 
the derivative terms on the right hand side of (8)-(10) is replaced 
by a certain difference, then all derivative terms on the right hand 
side of (8)-(10) are replaced by the same difference operator. 

(e) The discretization must be exact, i.e. produce exactly the same 
result as the analogous continuous operation, on all at least arbi- 
trary linear functions, i.e. f(x) = ax + b, a, b, arbitrary. 
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With these restrictions, it is fairly simple to show that the 
quasi-discrete schemes we wish to consider can all be reduced to the 



general form: 



oT ^ L cJ U ^ X,t ^ = f L g[v(x,t)] - g L^[h(x,t) ] 

ft * L c J V ^ X ’ t - ) ^ = ~ f [u(x , t) ] 



( 22 ) 



— (L [h(x , t) ] } = -H L [u(x , t) ] 
d t a a 



where [ ], L^[ ], and [ ] , are invariant linear discretization operators 

having the following general form: 

a 

L [f(x)] = a^f(x) 4- £ a^tfCx + nAx) + f(x - nAx) ] 



n=l 

n 



L [f(x)] = g Q f(x) + J B n [f(x + nAx) + f(x - nAx)] 



(23) 



n 



n=l 



[ f (x) ] = -^ l a n [f(x + nAx) - f(x - nAx)] 
n=l 



The various values of a , 3 , and a will be referred to as the filter 

n n n 

weights. Thus, for example, the finite difference mesh which Arakawa and 
Lamb (1977) call Scheme C (Figure 5) , and which is given in finite difference 
form as 



3u _ v(x+d/2,t) 4- v(x-d/2,t) _ h(x+d/2,t) - h(x-d/2,t) 

at " r 2 8 d 

av _ _ u(x+d/2,t) +v(x-d/2,t) 

at 2 

3h _ ,, u(x+d/2,t) - u(x-d/2,t) 

IF " " H a 

corresponds to the following filter weights: 

a Q = 1, B 1 = -|, = -| , all others zero. (Note that in this 

formulation Ax = d/2 .) 

Determination of the filter properties of the system represented by (22) 
requires taking the Fourier transform (in x) of those equations. However, the 
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shift property of these transforms (equation (9)) simplifies the terms 



involving the discretization operators, since, using (9) 



(LlfWl) -f I' c e ikn(ta >|f<k) • 



n=-n. 



n 



(24) 



Thus, the Fourier transform of (22) is 

. du 



ct(k)^- = f 0(k)v - i g a(k)h , 



a(k)^- = -f 3(k)u , 



a (k)^jjr = -i Ha (k)u , 



| 

t 



(25) 



where 



n 



a(k) = a + 2 I a cos(nkAx) , 
o n 

n=l 



(26) 



n 



e 

3 (k) = 3 Q + 2 l cos(nkAx) , 
n=l 



(27) 



n 

2 ^ 

a(k) = - — J a sin(nkAx) • 
Ax ^ n 
n=l 



(28) 



Equations (25) can be solved, using the same eigenvalue/eigenvector approach 
as was used for (12), to yield: 



a f 3v 

~ ~ A O 

u = u cos vt H sin vt - 

o v v 



A iCJ 8 h O A 

sin vt 



f gu 



v = - 



o . A 
sin v 



t + cos vt}v Q + (1 - cos vt}h Q 

v v v 



iaHu 



h = -- 



o . a a iaHfgr*. 

~ sin vt {1 - cos vtjv + 



I fV , a 2 gH 



cos vt/6 



where now 



a2 f 2 6 2 (k) + gHa 2 (k) 
v - 2 , 

a (k) 
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(29) 



or 



A 

V 



a(k) 



/? 



(k) 



2 2 

+ A a (k), 



v(k) 

o(k) ’ 



(30) 



where 

v(k) = f/g 2 (k) + X 2 a 2 (k) . (31) 

Comparison of (29)-(31) with ( 14 ) — ( 15) show that the main impacts 
of the quasi-discretization described by (22) are: 

(1) The filter coefficients 

1 k _k 

v ’ v ’ 2 

v 

and the squares of the first two, are replaced, respectively, 

g o gg 

v 5 v 5 2 

v 

and the squares of the first two of these expressions. 

(2) The phase velocity 

* = L / 22 

k k A + rk z 



is replaced by 



ka (k) 



A 2 



g + 



2 2 
X o 



Since the filter coefficients discussed above determine how the 
energy in the initial disturbance is partitioned into the transient and steady 
state fields, then the degree to which the discrete filter coefficients approxi- 
mate the continuous ones can be viewed as a measure of the degree to which the 
discretization accurately protrays the continuous model energy distribution. 

Any difference between a continuous (differential) coefficient, and the corres- 
ponding discrete coefficient, should be considered as introducing a distor- 
tion between the continuous and discrete solutions. 
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Similarly, any differences between the continuous phase and group velocities 
and the discrete ones will result in distortion due to spreading of the waves. 
This effect will be especially noticeable if the distortion is pronounced at 
(spatial) frequencies that carry significant energy, i.e. for which filter 
coefficients are "large". 

Actually, as noted before, the phase velocity is less important than 
the group velocity, 

dv 

dk 



in the continuous case, and 




in the discrete case, since this quantity measures the velocity at which the 
energy in the different frequencies propagates. 
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VI. COMPARISON OF THE ARAKAWA SCHEMES 



The specific discretization schemes which we shall analyze using 
(22)-(30) are based on the arrangements of sample points described in 
Arakawa and Lamb (1977) as Schemes A-D. The arrangements of these points 
in one dimension is shown is Figure 5. We shall consider both finite 
difference and finite element formulations of the discrete problem, 
whereas Arakawa and Lamb considered only second order finite differences. 

The finite difference formulations are derived for both second and fourth 
order approximations to the spatial derivative terms, while the finite 
element formulations are derived (for Schemes A-C only) using piecewise 
linear elements (the so-called "hat" function basis). The filter weights 
derived for each of these methods are shown in Tables 1-3. Note that since 
each scheme "samples" the initial disturbances at a distance of d, then the 
shortest wave that can be resolved in each case is of frequency (tt / d) . 

In Figures 6-8, the amplitude coefficients of the discretized approxi- 
mate "filters" are compared to the "filters" in the differential model, for 
the values of the physical parameters given in Arakawa and Lamb (1977). The 
following qualitative generalizations are clear: 

(1) In each of the methods (second order difference, fourth order 
difference, or finite elements) , the arrangement of points specified by 
Scheme A appears to produce the greatest distortion, and, furthermore, this 
distortion is most pronounced at high frequencies (short waves). 

(2) In each of the methods, Schemes B, C and D very accurately approx- 
imate the transfer function for the "high pass" filter. 

(3) In each of the methods. Scheme C understates the amount of energy 
distributed into the short waves by both the low-pass and band pass filters, 
while Scheme B overstates this. 
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(4) In all cases, the finite element formulation appears more 
"accurate 11 than fourth order differences, which in turn is more "accurate" 
than the second order differences. 

(5) The (group) velocity distortion is especially serious in Schemes 

A and D, since for short waves there is actually a reversal in the direction 
of propagation. 

Note that even though Scheme B seems to have the least "distortion" in 
each case, there is still significant phase and group velocity distortion 
near the cut-off. Thus, we should not conclude immediately that Scheme B 
will produce the least "noisy" solutions. This is because the apparent noise 
in the solution is a combination of both amplitude and phase behavior, and 
any phase distortion may lead to perceived "noise" in the solution if (as 
is the case in Scheme B) significant energy is retained (due to the amplitude 
distortion effects) in the short waves. In fact, it may be preferable to 
tolerate more amplitude distortion, if it acts to reduce the amount of 
energy in the short waves, which are generally the slowest propagating, and 
which produce the appearance of noise in the solutions. 

As a final check on our analysis, we compared the height fields for 
both the continuous and all the discrete models, using the same methodology 
as in Arakawa and Lamb (1977). An initial disturbance given by: 



u_ L (x,0) = 0 



h^x.O) = 0 



v i (x,0) = 



Vq , - a < x < a 
0 , otherwise 



was considered, with the physical parameters having the following values 

g = 10 m sec 



H = 10 m 
f = 10 -4 sec" 1 



a = d = X/2. 

- 18 - 



The transform of the initial v field 

^ 2Vq sin(ka) 



v. = 
1 



7T TT 

was then band limited to the region - -j- < k < — , and the exact solution 



for h was computed as the inverse transform of 

ikHf 



2V sin(ka) 

{1 - cos vt}{ : } 



in the continuous (differential) case, and the inverse transform of 

io B fH „ ,V _ % , r 2V 0 sin(ka \ 

- {1 - cos( ? t)H k } 



V 

in each of the discrete cases. In each instance, the inversion was accom- 
plished using Simpson’s rule on the interval (0, — ) , with 600 subintervals. 

The different h fields are shown at Figure 9. Note that these agree 
with our predictions, based on the filter analysis, in that: 

(1) In the second order finite difference methods. Scheme A displays 

a great deal of numerical noise, as predicted due to the great amplitude dis- 
tortion of its filter near (tt /d ) . Schemes B and D show somewhat less noise 

than Scheme A, but still more than Scheme C, this due to the manner in which 

2 

Scheme C controls the amplitude near (tt /d) in both the (a/V) and (ga/V ) terms 

(2) Similar comments hold for the fourth order finite difference methods 
though, in general, due to better phase propagation, the spreading of the 
energy in the short waves is not so pronounced, and hence the noise in 
Schemes B and D is reduced. 

(3) In the finite element methods, Scheme A still remains quite noisy, 
due to its amplitude distortion near (tt /d ) . Schemes B and C are virtually 
identical, since Scheme B seems to compensate for slightly more energy in 
the short waves by maintaining better phase relationship. 
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VII CONCLUSIONS 



In this paper we have adapted concepts that are commonly associated 
with the design of filters in electrical engineering to the analysis of 
discretized representations of partial differential equations, and spec- 
ifically to the one dimensional shallow water equations. We have shown 
that these concepts can provide valuable insights into why certain solution 
schemes provide "better" solutions than other, and a useful tool for com- 
paring both qualitative and quantative aspects of different schemes. 

We have clearly demonstrated that, in contrast to the analysis used in 
earlier studies, a complete analysis of the suitability of a certain discrete 
scheme to model a physical process, in our case the process of geostrophic 
adjustment, depends not only on providing reasonably accurate duplication of 
the phase propagation (delay) characteristics of the modeled system, but 
also on producing reasonably accurate duplication of the amplitude response 
characteristics. Furthermore, because of the phase distortion near the 
Nyquist cut-off of discrete schemes, it is often desirable to control the 
amplitudes of short waves. 

In the process of geostrophic adjustment studies, we have concluded that 
schemes that use unstaggered grid points are generally poor, irrespective of 
the particular finite difference/finite element method used, due to some in- 
herent properties of the model. The difficulties which have been reported 
with these methods in the literature, and especially the problem of "noise" 
in finite element models, appears directly attributable to the amplitude 
response near the Nyquist cut-off methods based on an unstaggered arrange- 
ment of points. This is certainly related to the tendency of schemes based 
on unstaggered grid points to produce solution separation, and appears to 
arise from the fact that the coupling of the fields in the basic differential 
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equation simply makes specifying all the quantities at every grid point an 
overspecification of the problem. 

Lastly, we conclude, that, based on the one dimensional evidence only, 
the method identified by Arakawa and Lamb (1977) as Scheme C has been quite 
popular, not only due to its good phase propagation characteristics, but 
also due to its tendency, in either finite difference or finite element for- 
mulations, to "window" out much of the high frequency noise in the discre- 
tized model. 

In summary, the concepts we have applied yield extremely valuable 
insights into the qualitative behavior of discrete methods for approximating 
partial differential equations. In many ways, we have only started to bring 
much of the emerging understanding of digital filters to bear on this problem. 
A significant, and interesting outgrowth will be the consideration of two 
dimensional problems. 
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TABLE 1 



FILTER WEIGHTS FOR SECOND ORDER FINITE DIFFERENCES 

(Ax = d/2) 



SCHEME 



DENOTED 



A2 



B 



B2 



1 



1 



0 



1 

2 



C 



C2 



1 



0 



1 1 
2 2 



D D2 



1 0 



1 

2 



0 
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TABLE 2 



FILTER WEIGHTS FOR FOURTH ORDER FINITE DIFFERENCES 

(Ax = d/2) 



SCHEME 



DENOTED 



A 



A4 



110 0 




J 

24 



B 



B4 



1 1 



0 



27 

48 



0 



_1 

48 



0 



C 



C4 



1 0 



1 _27 

2 48 



0 



JL 

48 



0 



D 



D4 



1 0 



1 

2 



0 




JL 

24 
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FILTER WEIGHTS FOR PIECEWISE LINEAR FINITE ELEMENTS 
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